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Abstract 

Monte Carlo simulations of lattice QCD at non-zero baryon chemical 
potential fi suffer from the notorious complex action problem. We con- 
sider QCD with static quarks coupled to a large chemical potential. This 
leaves us with an SU(3) Yang-Mills theory with a complex action con- 
taining the Polyakov loop. Close to the deconfmement phase transition 
the qualitative features of this theory, in particular its Z(3) symmetry 
properties, are captured by the 3-d 3-state Potts model. We solve the 
complex action problem in the Potts model by using a cluster algorithm. 
The improved estimator for the //-dependent part of the Boltzmann fac- 
tor is real and positive and is used for importance sampling. We localize 
the critical endpoint of the first order deconfmement phase transition line 
and find consistency with universal 3-d Ising behavior. We also calcu- 
late the static quark-quark, quark-anti-quark, and anti-quark-anti-quark 
potentials which show screening as expected for a system with non-zero 
baryon density. 
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1 Introduction and Summary 



Non-perturbative dense QCD can presently not be studied from first principles be- 
cause Monte Carlo simulations of lattice QCD with non-zero baryon chemical po- 
tential \i suffer from a severe complex action problem. The Boltzmann factor in the 
path integral can then not be interpreted as a probability and standard importance 
sampling methods fail. In particular, when the /i-dependent part of the Boltzmann 
factor is included in the measured observables, due to severe cancellations the re- 
quired statistics is exponentially large in the space-time volume |2|. 

The complex action problem prevents the numerical simulation of a large class of 
interesting physical systems including other field theories at non-zero chemical po- 
tential or non-zero ^-vacuum angle as well as some fermionic field theories with an 
odd number of flavors. A special case of the complex action problem is the so-called 
fermion sign problem which arises for fermionic path integrals formulated in a Fock 
state basis. The problem is due to paths that correspond to an odd permutation of 
fermion positions which contribute negatively to the path integral. There are nu- 
merous condensed matter systems ranging from the repulsive Hubbard model away 
from half-filling to antiferromagnetic quantum spin systems in an external magnetic 
field that cannot be simulated with standard Monte Carlo algorithms. Meron-cluster 
algorithms have been used to solve the sign or complex action problems in several 
of these cases. For example, the first meron-cluster algorithm has led to a solution 
of the complex action problem in the 2-d 0(3) symmetric field theory at non-zero 
vacuum angle ||. In this model, some of the clusters are half-instantons, so they 
are called meron-clusters. The complex action problem also arises in the 2-d 0(3) 
model at non-zero chemical potential. When formulated as a D-theory [f|, [5|, ^ [7j - 
i.e. in terms of discrete variables that undergo dimensional reduction — the complex 
action problem has also been solved with a meron-cluster algorithm |J. Recently, 
the meron concept has been generalized to fermions ||. Meron-cluster algorithms 
have led to a complete solution of the fermion sign problem in a variety of mod- 
els including non-relativistic spinless fermions [LTJ, relativistic staggered fermions 
n, 12, and some models in the Hubbard model family |], [14]].[] Recently, a 



meron-cluster algorithm has been used to solve the sign problem that arises for 



quantum antiferromagnets in an external magnetic field [15]. For a review of these 



recent developments and a preliminary version of the present results see [ lq , |17 . 



In the conventional formulation of lattice QCD the quarks are represented by 
Grassmann fields. When the quarks are integrated out, they leave behind a fermion 
determinant that acts as a non-local effective action for the gluons. At zero chemi- 
cal potential and for an even number of flavors, the fermion determinant is real and 
positive and can thus be interpreted as a probability for generating gluon field con- 
figurations. Despite the fact that standard importance sampling techniques apply, 



1 The models investigated so far only show s-wave superconductivity. 



2 



the non-local nature of the effective gluon action makes lattice QCD simulations 
with dynamical fermions very time consuming. With a non-zero chemical potential 
for the baryon number, the fermion determinant becomes complex and standard 
importance sampling techniques fail completely |2|. This is the reason why non- 
perturbative QCD at non-zero baryon density can presently not be studied from 
first principles. 

It is natural to ask if a meron-cluster algorithm could be used to solve the complex 
action problem in QCD. When one integrates out light quarks, one obtains a non- 
local effective action for the gluons and it appears unlikely that the meron concept 
will apply. On the other hand, when one describes the quarks in a Fock state basis, 
the complex action problem is still present, in the form of a fermion sign problem. 
Our hope is that this problem will eventually be solved by a meron-cluster algorithm 
applied to the D-theory formulation of QCD @, [5|, |], [7| , since the quark and gluon 
degrees of freedom are then discrete and are much easier to handle. In this paper 
we address a simpler situation first. We consider QCD in the limit of very heavy 
quarks with a large chemical potential. These can be integrated out, introducing 
Polyakov loops into the effective gluon action. When quarks are integrated out at 
non-zero chemical potential fi we expect a complex action, and in this case it arises 
because a Polyakov loop $ and its charge conjugate $* get different weights when 

Polyakov loops are only non-local in the Euclidean time direction, so this effective 
gluon action is more manageable than the one that arises for a general fermion 
determinant. Indeed, Blum, Hetrick and Toussaint have simulated the theory in 
this form on lattices of moderate size where the complex action problem is less 
severe Recently, Engels, Kaczmarek, Karsch and Laermann have studied QCD 
with heavy quarks at fixed baryon number. Again, for moderate baryon density and 
moderate volumes the complex action problem is not too severe and simulations 
are possible |19|| . Ultimately one would like to be able to solve the complex action 



problem for this gluon action completely. At the moment, we still cannot apply a 
meron-cluster algorithm to solve the problem, because the construction of efficient 
cluster algorithms for non-Abelian gauge theories seems to be impossible for Wilson's 
formulation of lattice field theory. Here we will simplify the problem further by 
replacing the gauge dynamics by that of the Z(3) Potts model representing the 
Polyakov loops p0, |21| . We have found a cluster algorithm that solves this complex 



action problem in the Potts model approximation to QCD. 

The 3-d Z(3)-symmetric Potts model has often been used as an approximation to 
QCD with static quarks. In particular, the phase transition to a broken Z(3) sym- 
metry phase at high temperature corresponds to the first order deconfining phase 
transition in QCD. As has been noted by Condella and DeTar, a term that corre- 
sponds to a chemical potential can also be included in the Potts model, explicitly 
breaking the Z(3) symmetry ||22|| . As the coefficient of this term grows, the first 
order deconfinement phase transition persists but it becomes weaker and ultimately 
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Figure 1: The phase diagram of the Z(3) Potts model in the (h, k) -plane. The 
ordinary deconfinement phase transition at T = (0,0.550565(10)) is a triple point 
from which a line of first order phase transitions emerges. This line terminates 
in the critical endpoint E = (0.000470(2), 0.549463(13)) and continues only as a 
crossover. 

disappears in a critical endpoint. This point is expected to be in the universality 
class of the 3-d Ising model. In this paper we will confirm this expectation with 
numerical simulations. 

In principle, one can imagine deriving an effective 3-d 3-state Potts model directly 
from QCD by integrating out all degrees of freedom except for the Z(3) phase of the 
Polyakov loop. However, the resulting Potts model action would be very complicated 
and cannot be derived in practice, except in the strong coupling limit. Here we 
approximate QCD with heavy quarks by a 3-d Z(3)-symmetric Potts model with 
a standard nearest-neighbor interaction. Universal features like the nature of the 
critical endpoint of the deconfinement phase transition are correctly reproduced in 
this approximation. Figure 1 contains the phase diagram of the 3-d 3-state Potts 
model in the (h, «)-plane. The parameter h represents exp (/?(// — M)) in QCD with 
quarks of mass M at chemical potential //. We study the limit M, /j, — > oo for any 
given fi — M. Large h corresponds to fi > M and small h to fi < M. Because 
fi — M <C M, ji we are always, for any h, in the immediate neighborhood of the 
onset of non-zero density for the heavy quarks. This means that it does not matter 
whether they are fermions or bosons, since they never move. The difference would 
only become apparent above the onset, where either a Fermi surface or a degenerate 
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Bose gas would occur, and our order of limits is such that we never get that far from 
the onset. The parameter k is the standard Potts model coupling, which corresponds 
roughly to the temperature T = 1/(3. The ordinary first-order deconfinement phase 
transition at h — (point T in Figure 1) extends into a line of first order transitions 
that terminates in the critical endpoint E. This endpoint occurs at such a low value 
of h that the complex action problem is not very severe there, and we found that 
the most efficient way to locate and study it was to employ a reweighted Metropolis 
algorithm, which can in this case be applied at volumes large enough to show the 
critical behavior. Similar methods were used recently by Karsch and Stickan |^3| in 
a version of the 3-d 3-state Potts model where the action is real, and the endpoint 
was found to have the critical exponents of the 3-d Ising model. We find that in 
the Potts model with complex action the endpoint has the same critical properties. 
Furthermore its position is barely shifted in comparison to the model with real 
action. In this paper we do not limit our attention to the endpoint, but develop a 
method that solves the complex action problem everywhere in the phase diagram. 

We also calculate the potentials between static quarks and anti-quarks in the 
Potts model approximation to QCD. In the confined phase at fi = the static 
quark- ant i- quark potential is linearly rising with the distance as a manifestation of 
confinement. For the same reason the quark-quark and anti-quark-anti-quark po- 
tentials are infinite at all distances. In the deconfined phase the quark-anti-quark 
potential reaches a plateau at twice the (now finite) free energy of a quark. Simi- 
larly, the quark-quark and anti-quark-anti-quark potentials are no longer infinite. It 
should be noted that quark-quark and anti-quark-anti-quark potentials are usually 
not calculated in lattice simulations. This is because — as a consequence of the Z(3) 
Gauss law — quark or anti-quark pairs cannot exist in a finite spatial volume with 
periodic boundary conditions [p4[ . Interestingly, this changes for ^ because 



then there are compensating background charges in the medium that can absorb 
the Z(3) flux of an external quark. Since the chemical potential explicitly breaks 
the Z(3) symmetry, there is no longer a clear distinction between confinement and 
deconfinement for /j ^ 0. This manifests itself in the phase diagram by the fact 
that confined and deconfined phases are analytically connected. Figure 2 shows the 
quark-anti-quark, quark-quark and anti-quark-anti-quark potentials on the confined 
side (a) and on the deconfined side (b) of the crossover. Note that at ji ^ even 
in the confined phase the quark-anti-quark potential now reaches a plateau. The 
plateau height corresponds to the sum of the free energies Fq of an external quark 
and Fq of an external anti-quark. For /i > quarks are favored in the medium 
while anti-quarks are suppressed. As a consequence, the free energy of an external 
static quark is larger than that of an external static anti-quark. While an external 
static anti-quark can bind with a single background quark from the medium and 
form a meson, an external static quark needs two quarks from the medium to form 
a baryon. Indeed, on the confined side of the transition Fq is clearly larger than Fq, 
while on the deconfined side Fq and Fq are more or less the same. We have nor- 
malized the potentials such that at zero distance a static quark-anti-quark pair has 
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Figure 2: The static quark- anti- quark, quark-quark and anti- quark- anti- quark poten- 
tials (a) on the confined side (at h = 0.01, k = 0.50) and (b) on the deconfined side 
(at h = 0.01 and k = 0.56) of the crossover. 
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zero energy. In the Potts model, two quarks at zero distance are indistinguishable 
from a single anti-quark, and similarly, two anti-quarks on top of each other behave 
like a single quark. Hence, at zero distance the quark-quark potential Vqq(O) agrees 
with the free energy of a single anti-quark Fq and the anti-quark-anti-quark poten- 
tial obeys Vqq(0) = Fq. At asymptotic distances the potentials Vqq(oo), Vqq(oo) 
and Vqq(oo) take the values Fq + Fq, 2Fq and 2Fq, respectively. This behavior is 
consistent with our numerical data shown in fig. 2. 

In the absence of a chemical potential, the Potts model can be simulated with 
the original Swendsen-Wang cluster algorithm E5|. When a chemical potential is 
introduced, the Potts model suffers from the complex action problem and standard 
importance sampling methods including the cluster algorithm fail. In this paper, we 
construct an improved estimator for the /^-dependent part of the Boltzmann factor 
by averaging analytically over all configurations related to each other by cluster 
flips. In contrast to the original Boltzmann factor, the improved estimator is real 
and positive and can be used for importance sampling. This solves the complex 
action problem completely. 

Although the Potts model inherits the complex action problem from QCD, it 
can be transformed into a "flux model" that has no complex action problem fl22"f . 



The flux model has been simulated in p2| and the disappearance of the first order 



deconfining phase transition at large chemical potential has been observed numeri- 
cally. These results may give encouragement to the hope that QCD itself could be 
transformed into a model without a complex action problem. In this paper we show 
that, at least for the Potts model, it is more efficient to leave it in its usual form 
and solve the complex action problem with our cluster algorithm than to transform 
it into a flux model and use conventional Metropolis methods. 

The paper is organized as follows. Section 2 contains a derivation of the effective 
gluon action resulting from static quarks with large chemical potential as well as 
its Potts model approximation. In section 3 we describe the cluster algorithm that 
solves the complex action problem. Section 4 contains the derivation of the flux 
representation of the Potts model and a description of a Metropolis algorithm to 
simulate it. A comparison of the Metropolis algorithm for the flux model and the 
cluster algorithm for the original Potts model shows that the latter is more efficient. 
In section 5 we present the physical results concerning the critical endpoint E. Using 
finite-size scaling techniques, we are able to determine the position of the critical 
endpoint of the deconfinement phase transition to high accuracy. Our results are 
consistent with the expected universal 3-d Ising behavior. Finally, section 6 contains 
our conclusions. 
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2 QCD with Heavy Quarks and the 3-d 3-State 
Potts Model 



The partition function for a pure SU(3) Yang-Mills theory is given by 

Z = j VA exp(-S[A]), (2.1) 

where ^ 

S[A] = J dtj d 3 x J-^Tr[F^F^\, (2.2) 

is the Euclidean action for the gluons and (3 is the inverse temperature. The action 
is invariant under gauge transformations 

g(x,0) = g{x,P)z, (2.3) 

that are periodic in Euclidean time up to an element z of the center Z(3) = 
{exp(27un/3), n = 1,2,3} of the non-Abelian gauge group. In the presence of a 
single external heavy quark of bare mass M at an undetermined position x the 
partition function turns into 



Z Q = JvA$[A] exp(S[A])exp(-pM), 



(2.4) 



where 

$[A] = J d 3 x Tr[Pexp(-^ dt A 4 {x,t))}, (2.5) 

is the spatial integral of the Polyakov loop. Ultimately, the mass M will be sent to 
infinity. Note that while the center transformation of eq. (|2.3|) leaves the pure gluon 



action S[ 9 A] = S[A] invariant, the Polyakov loop transforms into 

$[ 9 A]=z$[A]. (2.6) 

This shows that in the presence of the external quark, the Z(3) symmetry is explicitly 
broken. The partition function for a system of gluons in the presence of a single 
heavy anti-quark is given by 

Z Q = JvA $L4]* exp(-S[A]) exp(-/?M), (2.7) 

where * denotes complex conjugation. Let us now consider a system of gluons in a 
background of n static quarks and n static anti-quarks. The partition function then 
takes the form 

Z n , n = [VA ^[A] n ^mA]*fex V (-S[A])ex V (-/3M(n + n)). (2.8) 
j n. n. 
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The factors 1/n! and 1/n! appear because quarks are indistinguishable, as are an- 
tiquarks. Introducing the quark chemical potential \i that couples to (n — n), i.e. 
three times the baryon number, we obtain the grand canonical partition function 

= [vA^^A] n ^($[AY) n exp(-S[A]-(3n(M-ii)-(3n(M + fi)) 
j n. n. 

n,n 

= J VA exp(-S[A] + exp(-p(M - n))$[A] + exp(-/3(M + /i))$[A]*). 

(2.9) 

As expected, the presence of quarks and anti-quarks leads to an explicit breaking of 
the Z(3) center symmetry. Furthermore, in the presence of a non-zero chemical po- 
tential the effective action for the gluons is complex. Note that in the £77(2) case the 
action remains real because then the Polyakov loop itself is real, i.e. $[v4]* = 3>[A]. 
The action becomes real even in the SU(3) case if \x is purely imaginary. Further- 
more, one can see that the chemical potential explicitly breaks the charge conju- 
gation symmetry that replaces <&[A] by $L4]*. In fact, under charge conjugation 
the action turns into its complex conjugate. We have assumed that the quarks are 
static. Hence, to be consistent we must consider the limit M — > oo. In order to 
obtain a non-trivial result, we simultaneously take the limit \x — > oo such that M — fi 
remains finite. The partition function then simplifies to 

Z(fj.) = j VA exp(-S{A] + exp(-/3(M - n))$[A]). (2.10) 

As discussed in [IS| and a similar result can be obtained by simplifying the full 
QCD quark determinant in the static quark limit. In general the determinant would 
contain all Wilson loops, but because M is large most of them are suppressed. The 
only ones that survive are those for which the enhancement due to the chemical 
potential compensates for the suppression due to the mass. These are the Polyakov 
loops that progress in a straight line from Euclidean time t — to t — {3 at some 
position x. In the loop expansion of the quark determinant, each of these has a 
weight exp(— (3{M — //)). 

Up to this point we have treated QCD consistently in the static quark limit. 
The resulting effective action for the gluons is complex and we presently don't know 
how to simulate it efficiently. For that reason we now replace the gluon system by 
a simple 3-d lattice 3-state Potts model. The Potts spins $ x e Z(3) replace the 
original Polyakov loop variables and the partition function turns into 



Z{h) = J X?$ exp(-S[$] + hJ2^x), (2.11) 



where h replaces exp(— (3{M — fx)). Note that the Potts model action is still complex. 
In principle, one can imagine integrating out all QCD degrees of freedom except for 
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the Z(3) phase of the Polyakov loop and thus derive an effective Potts model action 
directly from QCD. In practice this is impossible, except in the strong coupling 
limit. For simplicity, we therefore replace the pure gluon action S[A] by a standard 
nearest-neighbor Potts model interaction 



The coupling constant k is not related in a simple way to the parameters of QCD. 
Still, a large value of k corresponds qualitatively to the high-temperature deconfmed 
phase, while small k values correspond to the confined phase. As mentioned in the 
introduction, the Potts model also retains the general features of the QCD phase 
diagram. At h — (M infinite, \i finite) there is a first-order phase transition as 
a function of k, between the disordered (confined) phase that respects the Z(3) 
symmetry and the ordered (deconfmed) phase that spontaneously breaks it. An 
order parameter for this transition is ($) . As h rises from zero, the chemical potential 
term explicitly breaks the Z(3) symmetry, the phase transition weakens, and then 
ends at a critical point. Correspondingly, in heavy-quark QCD the quarks begin to 
contribute to the partition function when /i gets close to M, and there is no longer 
an order parameter for deconfinement. The deconfining phase transition terminates 
at a critical endpoint. 

3 Cluster Algorithm Solution of the Complex Ac- 
tion Problem 

In this section we first discuss the general nature of the complex action problem and 
then discuss the cluster algorithm that solves this problem for the Potts model. We 
also construct improved estimators for various physical quantities. 

3.1 The General Nature of the Complex Action Problem 

When the action is complex the resulting Boltzmann factor cannot be interpreted 
as a probability and hence standard importance sampling techniques fail. When one 
uses just the absolute value of the Boltzmann factor for importance sampling and 
includes its complex phase in measured observables O, expectation values take the 
form 




(2.12) 



X.l 




X 



(Oexp( y ihJ2J™$x))R 
{ex.p(ihJ2 x lm ®x))n 



(3.1) 
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The subscript R refers to 
partition function 

Z R 

By definition we have 
(exy{ihyirn<b x )) R 

X 

where / and f R are the free energy densities of the original complex and the modified 
real action systems, respectively, and V is the spatial volume. Hence, the denomina- 
tor in eq. (|3.1|) becomes exponentially small as one increases the volume. The same 
is true for the numerator, because (O) itself is not exponentially large in V. 

Although, in principle, simulating the modified ensemble is correct, in practice 
this method fails for large volumes. The reason is that observables are obtained 
as ratios of exponentially small numerators and denominators which are themselves 
averages of quantities of order one. This leads to very severe cancellations and 
requires an exponentially large statistics in order to obtain accurate results. To see 
this, we estimate the relative statistical error in the determination of the average 
phase of the Boltzmann factor exp(ih J2 X Ini^a,). Since (exp(ih J2 X ^- m ^x))R = Z/Z R 
the average itself is real. When one generates N statistically independent field 
configurations in a Monte Carlo simulation, the resulting error to signal ratio is 
given by 

Aexp(z/^,Im$ x ) _ ^/(\exp(ihJ2 x Im$ x ) - (exp(ihJ2 x Im$ x )) R \ 2 ) R 
{exp{ih^ x Im<S> x )) R VN(exp(ihJ2 x lm$ x )) R 

y/j - (exp(^E,Im$,))| ^ expjVtf - f R )) 
^(exp(ihJ2 x hn^ x )) R ~ VN 

For large V we have used (exp(ih J2 X bm^))/? C 1 as implied by eq.([3.3|). Conse- 
quently, in order to obtain an acceptable error to signal ratio one must generate at 
least N « exp(2V(f — f R )) configurations. For large volumes this is impossible in 
practice. 



a modified ensemble with a real action described by the 
= y X?$ exp(-£[$] + /i^Re$). (3.2) 



— / £>$ exp(i/i^Im$ :E )exp(-S'[$] + /iJ^Re$) 

R J X X 

A « exp(-V(/ - f R )), (3.3) 

^R 



3.2 The Cluster Algorithm for the Potts Model 

Let us now outline the ideas that underlie the cluster algorithm that we use to solve 
the complex action problem. It is based on the original Swendsen-Wang cluster 
algorithm [^5j for the Potts model without chemical potential. In fact, in the limit 
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h — our algorithm reduces to that algorithm. The Swendsen-Wang cluster algo- 
rithm decomposes the lattice into independent clusters of connected sites. Each spin 
belongs to exactly one cluster and all spins within a cluster are assigned the same 
random Z(3) element. In this paper, we construct an improved estimator for the 
/i-dependent part exp(h J2 x $ x ) of the Boltzmann factor by analytically averaging 
it over all configurations related to each other by cluster flips. Although, for an 
individual configuration exp(/i J2 X ®x) is m general complex, its improved estimator 
is always real and positive and can thus be used for importance sampling. This 
completely solves the complex action problem. 

Let us first describe the original Swendsen-Wang algorithm for h — 0. In this 
method one introduces variables b = 0, 1 for each bond connecting neighboring 
lattice sites x and y = x + i and one writes the nearest neighbor Boltzmann factor 

as 

exp(/c<J* B) * tf ) = ( e * - 1) + 5 b:0 \. (3.5) 

6=0,1 

In the enlarged configuration space of spin and bond variables, the bond variables 
impose constraints between the spin variables. When a bond is put (i.e. when 6=1), 
the spin Boltzmann factor is S& x ,& v (e K — 1) and hence the spin variables and $ y 
at the two ends of the bond must be identical. On the other hand, when the bond 
is not put (6 = 0), the spin Boltzmann factor is 1 and thus the variables & x and <3> y 
fluctuate independently. The spin variables, in turn, determine the probability to 
put a bond. When the spins $ x and Q y are different, the bond Boltzmann factor is 
5b t o an d thus the bond is not put. On the other hand, when <!> x and <£> y are the same, 
the bond Boltzmann factor is [5b,i(e K — 1) + 5j, j0 ]. Consequently, a bond between 
parallel spins is put with probability p—1 — e~ K . Note that for k — no bonds are 
put, while for k = oo parallel spins are always connected by a bond. 

The Swendsen-Wang cluster algorithm updates bond and spin variables in alter- 
nating order. First, for a given spin configuration, bonds are put with probability 
p between parallel neighboring spins. No bonds are put between non-parallel spins. 
Then the spins are updated according to the constraints represented by the resulting 
bond configuration. Spins connected by bonds must remain parallel, while spins not 
connected by bonds fluctuate independently. Hence, to update the spins, one must 
identify clusters, i.e. sets of spins that are connected by bonds. All spins in a cluster 
are parallel and are assigned the same random Z(3) element in the spin update. All 
spins belong to exactly one cluster. It should be noted that a cluster may consist of 
a single spin. A configuration consisting of Nq clusters can be viewed as a member 
of a sub-ensemble of 3^° equally probable configurations which result by assigning 
Z(3) elements to the various clusters in all possible ways. As was already pointed out 
by Swendsen and Wang, one can construct improved estimators for various physical 
quantities by averaging analytically over all 3^° configurations in a sub-ensemble. 
Since the number of clusters is proportional to the volume, this effectively increases 
the statistics by a factor that is exponentially large in V . 
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Let us construct an improved estimator for the /i-dependent part exp(/i^ a , <& x ) 
of the Boltzmann factor. Although for an individual configuration this term is in 
general complex, its average over a sub-ensemble of 3^° configurations is always real 
and positive. This results from the following observations. The /i-dependent part of 
the Boltzmann factor is a product of cluster contributions 

exp(hJ2^) = l[exp(hJ2^)- (3-6) 

x c iec 

Since the clusters are independent, the sub-ensemble average is a product 

<exp(/i 5^)3*0 = H(exp(hJ2®x)) 3 i (3.7) 

x C x<=C 

of 3-state averages for the individual clusters 



3 

x&c $ez(3) 

= ^[exp(/i|C|) + 2exp(-/i|C|/2)cos(v / 3/i|C|/2)] 

= W(C), (3.8) 

which defines a weight W(C) for each cluster. We have used the fact that all spins 
$ x in a given cluster C take the same value $ G Z(3) so that YlxeC ® x = 1^1 ^ 
where \C\ = XLec ^ * s ^ ne cms t er si ze - It is easy to show that the expression in 
eq.( |3.8| ) is always positive and can hence be used for importance sampling. This is 
crucial for a complete solution of the complex action problem. 

For a given bond configuration one can integrate out the spin variables and one 
obtains 

Z = I Vb(e K - l) Nb 3 Nc Y[ W(C) 
J c 

Vb (e K - 1)^ JJ[exp(/i|C|) + 2exp(-/i|C|/2) cos(V3h\C\/2)]. (3.9) 



o 



Here Nb is the number of bonds that are put (i.e. have 6=1). The factor 3^° 
represents the number of allowed spin configurations for a given bond configuration 
and the factors W(C) come from the improved estimator. The effective action for 
the bond variables depends only on the sizes |C| of the clusters corresponding to a 
given bond configuration. Note that the factor 1/3 per cluster in eq. (|3.8| ) cancels 
against the factor 3 Nc . 



Our algorithm directly updates the partition function of eq . ( p79|) , i.e. it only op- 
erates on the bond variables while the spins are already integrated out analytically.^] 

2 B. Scarlet was first to realize that the spin variables need not even be simulated. 
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The bond variables that define the clusters are updated with a local algorithm. A 
bond whose value does not affect the cluster sizes is put with probability p = l — e~ K . 
This happens when the two sites at its ends belong to the same cluster because they 
are connected indirectly through other bonds. A bond whose value affects the clus- 
ter sizes is put with a probability that depends on the sizes of the corresponding 
clusters. When the bond is not put (6 = 0), its endpoints x and y belong to two 
different clusters C\ and C2 of sizes \C\\ and |C 2 | and the corresponding Boltzmann 
weight is 3 2 W (C\)W '(C2) • On the other hand, when the bond is put (b = 1), its 
endpoints belong to the combined cluster C\ UC2 of size \Ci \ + |C 2 | . In that case, the 
Boltzmann weight is 3W(Ci U C2){e K — 1). Hence, the bond is put with probability 

W(C^C^-l) 



3W , (Ci)W(C 2 ) + W(d U C 2 )(e K - 1) ' 



3.3 Improved Estimators for Physical Quantities 

In order to measure physical observables, it is crucial to construct improved estima- 
tors for them as well. Here we construct improved estimators for the Polyakov loop 
$ x , its charge conjugate $*, as well as for the correlators &x$y and $l&* y - 

The expectation values 

(cty = exp(-/3F Q ), ($:) = exp(-/3F ), (3.11) 

determine the free energies Fq of a quark and Fq of an anti-quark. The Polyakov 
loop correlators determine the quark- ant i- quark potential Vqq(x — y), the quark- 
quark potential Vqq(x — y) and the ant i- quark- ant i- quark potential Vqq(x — y) via 

exp(-f3V QQ (x-y)) = (^;), 
ex.p(-pV QQ (x - y)) = 

exp(-(3V QQ (x-y)) = (3.12) 
The improved estimator for the Polyakov loop is given by the sub-ensemble average 

(ct^exp^cty^ =i J2 $eMh\C x m J] W(C), (3.13) 

: *ez(3) c+c x 

where C x is the cluster that contains the point x. Hence, we obtain 

( ^ ) = ^/ P6 3WT ^ *e*p(h\C x \*)(e!<- 1) N > 3 N °l[W(C), (3.14) 

^ ' J ^ x > $ez(3) C 

i.e. after integrating out the spin variables, the Polyakov loop is represented by 

® X = W(C~) £ $ex P(^l $ )- (3.15) 

^ X > $GZ(3) 
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Similarly, the operator representing the charge conjugate Polyakov loop is 

V X > $GZ(3) 

The improved estimator for the Polyakov loop correlator & x &* is given by the sub- 
ensemble average 

(ct^exp^cty^ = I ®eMh\C x \$) ± E $*exp(/i|C7 w |$) 

: $GZ(3) $*eZ(3) 

x n w (°^ ( 3 - i? ) 

if the points a; and y belong to two different clusters C x and C y . If the points x 
and y belong to the same cluster (i.e. if C x = C y ), the improved estimator is simply 
given by 

($ x $; exp(/i J2 $ *)>3"c =]JW(C), (3.18) 

2 C* 

because then § X Q* = 1. Hence, the operator representing the correlator is 

= QWr L n E *exp(h\C x \*) E ^exp(fc|C tf |#), if C„ 
914/(014/(0,)^^ $gz(3) 

= 1, if C« = C y . (3.19) 
Similarly, we have 

= 9VTfC \\V(C ) S $ex P(^l $ ) E ^exp^C^), if C X ^C W , 

I V 2// $GZ(3) $eZ(3) 

= ^77?M E * 2 exp(/i|C x |<t>), ifC7 x = C7 y , 
^o*; $ ez(3) 

K% = QW (c\w(C) ^ $ * ex P(^l*) E S'expWC,,!*), if C^C S , 

^ av I yj $ G Z(3) $GZ(3) 

^ = E (**) 2 exp(/i|C x |ct>), ifC7 x = C7 w . (3.20) 

^ x > $ez(3) 



3.4 Severity of the Complex Action Problem 

In order to estimate the severity of the complex action problem, we like to determine 
the expectation value of the complex phase of the Boltzmann factor in the modified 
real action ensemble 

(exp^Vlm^))^^. (3.21) 
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Figure 3: The expectation value of the phase factor as a function of the volume at 
the critical endpoint E. 

Rather than implementing this directly in a simulation that uses the absolute value 
of the Boltzmann factor for importance sampling, one can measure Zr/Z with the 
cluster algorithm. In fact, an improved estimator for this quantity is given by 
YI C Wr(C)/W(C), where 

W R {C) = (exp^^Re^a = \[exj>(h\C\) + 2 exp(-h\C\/2)], (3.22) 
xec 6 

is the weight that replaces W(C) in the real action ensemble. Alternatively, one can 
construct a cluster algorithm that simulates the real action ensemble. In that case, 
one needs to measure Ylc W(C)/Wr(C) in order to obtain Z/Zr. 

Figure 3 shows (exp(ih ^ x Iui$> x ))r as a function of the volume V = L 3 at the 
critical endpoint of the transition line (point E in figure 1). Indeed, one finds an 
exponentially small signal, as expected from eq. (|3.3j ). Defining a scale parameter L 
that measures the severity of the complex action problem by 

(exp{ih J2 ^x))r oc exp ^-^Q (3.23) 

we get Lq pa 80 at the endpoint E. This means that the complex action problem 
at the endpoint is extremely mild. In fact, in practice it is not a problem at all 
up to volumes as big as 100 3 . Since the current computer hardware restricts the 
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Figure 4: The scale parameter Lq related to the severity of the complex action prob- 
lem as a function of h at k = 0.5495. The solid line is a spline to guide the eye. 
The phase transition takes place at the dotted vertical line. 

system size to a couple of million degrees of freedom, one can also study the point 
E with an algorithm that does not solve the complex action problem. On very large 
lattices the meron-cluster algorithm will become superior also at the point E since 
its computational effort is polynomial in the system size as opposed to exponential 
for the reweighted Metropolis algorithm. 

It should be noted that the complex action problem is most severe for interme- 
diate values of h. While it is obvious that there is no complex action problem at 
h — 0, it is perhaps less obvious that there is also no problem for large h. This is 
because 



approaches 1 in the limit h — > oo so that (exp(ih^2 x lm^ x )ji)ji = Zr/Z — > 1. Figure 
4 shows the complex action problem scale parameter Lq as a function of h for fixed 
k = 0.5495. It has a minimum at h ~ 1 meaning that the complex action problem 
ist most severe in that region. Defining the "practical complex action problem" as 
being present when {exp(ihJ2 x Im& x )) R < 0.01 on a 100 3 lattice (i.e. L < 60) we 
see that at k — 0.5495 there is no practical complex action problem for h < 0.003 
as well as for h > 6. A more physical definition of a practical complex action 
problem would compare L with the correlation length £ but we have not measured 




(3.24) 



W(C) exp(h\C\) + 2 exp(-h\C\/2) cos(v^/i|C|/2) 
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the correlation length. 



4 Flux Representation of the Potts Model 

In this section we map the Potts model to an equivalent flux model that does not 
suffer from the complex action problem. Then we describe a Metropolis algorithm 
to update the flux model and compare its efficiency with the cluster algorithm. 



4.1 Mapping the Potts Model to a Flux Model 



As pointed out in |E| the Potts model can be rewritten as a flux model that does 
not suffer from the complex action problem. As we have seen, the complex action 
problem can also be solved in the cluster formulation. Hence, the question arises 
if the flux or the cluster formulation leads to more efficient numerical simulations. 
Let us first match the flux model to the original Potts model. The flux model is 
formulated in terms of "electric" charges Q x G {0, ±1} defined on the lattice sites 
x and electric flux variables E x ^ e {0, ±1} living on the links. The charge and flux 
variables are related by the Z(3) Gauss law constraint 



Q x = J2(E x ,i - Earned 3. (4.1) 

i 

The action of the flux model takes the form 

S[E, Q] = £ E l + P E^ - ( 4 ' 2 ) 

x,i x 

The mass M and the chemical potential /i of the Z(3) charges are not directly related 
to the mass and chemical potential of quarks in QCD (which we also denoted by M 
and \i in section 2) but qualitatively they play the same role. The partition function 
of the flux model takes the form 

Z = U E II E U 5 ^M-S[E,Q]). (4.3) 

x Q x e{0,±l} x,i E Xii e{0,±l} x 

The 5-function 5 X imposes the Gauss law of eq. ([4.1|) at the point x and can be 
written as 

Inserting this as well as eq.( [4.2| ) for the action in eq.( [4.3|) one can integrate out the 
E x> i and Q x variables. The result of the E x i integration is 

E ^,.,)' : " exp (-y^) = 1 + 2exp I'^'i'I'./I*.,. . ,)• (4.5) 
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In the Potts model (up to an overall factor) the corresponding term is exp(K8<z> xi <s>^ ,). 
Thus, the flux model matches the Potts model if 

l + 2exp(-fl 2 /2) 

exp ( K ) = ; — — — . 4.6 

U ^ 1 - exp(-<772) 1 ; 

Note that the g — > limit of the flux model corresponds to the k — > oo limit of the 
Potts model. When one integrates out the charges Q x one obtains 

®f eM-MQ 2 x - /iQ x ) = 1 + exp(-/?(M - + exp(-/3(M + 

Q X £{0,±1} 

(4.7) 

In the original Potts model (up to an overall factor A) the corresponding term is 
exp(/i$ a; ). Hence, the flux model matches the Potts model if 



- exp(/i$) = A, 

$GZ(3) 

- ^2 $exp(/i$) = Aexp(-P(M + fj)), 



3 

*ez(3) 



- ®* ex P(^) — Aexp(—{3(M - /!)). (4. 



3 

$GZ(3) 



These relations can be used to determine the parameters exp(— (3(M — //)) and 
exp(— f3(M + /i)) of the flux model in terms of the Potts model parameter h. 



4.2 Metropolis Algorithm for the Flux Model and Compar- 
ison with the Cluster Algorithm 



As first described in ||22|| , the flux model can be updated with a simple Metropolis al- 
gorithm. One basic move in the algorithm creates or annihilates a nearest-neighbor 
charge- ant i- charge pair across a given link. The other basic move creates or annihi- 
lates an electric flux loop around an elementary plaquette. These basic moves are 
proposed on every link and plaquette and are accepted or rejected in a Metropo- 
lis step. We have implemented both the Metropolis algorithm for the flux model 
and the cluster algorithm for the Potts model and we have verified that physical 
observables obtained with the two algorithms agree with each other. 

The question arises which of the two algorithms is more efficient. The Metropolis 
algorithm is expected to suffer from critical slowing down at the endpoint of the first 
order phase transition with a dynamical critical exponent z ~ 2. Cluster algorithms 
are known to drastically reduce critical slowing down, in some cases even to z ~ 0. 
However, our cluster algorithm cannot eliminate critical slowing down completely 
because the decision to put a bond is more time-consuming than the one in the 
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Swendsen-Wang algorithm. For example, when a bond is removed, one must check 
if an old cluster decomposes into two new clusters. To minimize the computational 
effort, we simultaneously grow two clusters from the two ends of the bond. Once 
they touch each other, we know that the old cluster did not decay. Still, in the less 
likely event that the old cluster does decay, one must completely grow the smaller 
of the two clusters in order to decide if the bond can be deleted. Here we do not 
attempt to determine the dynamical critical exponent z of our cluster algorithm at 
the critical endpoint. 

An even more severe super-critical slowing down is expected close to the first 
order phase transition line. Along that line, the deconfined phase coexists with the 
confined phase, and the Monte Carlo simulation must tunnel between the two phases. 
In order to tunnel between the confined and the deconfined phase, a local algorithm 
must go through configurations containing both phases simultaneously. Since the 
interface that separates the two phases has non-zero interface tension, such configu- 
rations are exponentially suppressed. Hence, close to a first order phase transition, 
a local algorithm like the Metropolis algorithm necessarily suffers from exponential 
slowing down. This super-critical slowing down is even more severe than the power- 
law critical slowing down at a second order phase transition. Hence, we expect that 
the Metropolis algorithm for the flux model is not very well suited to study the first 
order phase transition line. In some cases, cluster algorithms can even eliminate 
super-critical slowing down. For example, in the broken phase of the Potts model 
three distinct deconfined phases coexist with each other. The Swendsen-Wang algo- 
rithm can efficiently tunnel from one deconfined phase to another because it assigns 
the same random Z(3) element to all spins in a cluster in a non-local spin update. 
Still, even the Swendsen-Wang algorithm suffers from super-critical slowing down 
at the first order phase transition that separates the confined from the deconfined 
phase. Although cluster flips can naturally lead to tunneling between distinct de- 
confined phases, they do not lead directly from a deconfined to the confined phase. 
To cure this problem, Rummukainen |26| has combined the cluster algorithm with 



the multi-canonical methods of Berg and Neuhaus |27] which can reduce the ex- 
ponential super-critical slowing down to a power-law behavior. Although this may 
well be possible, we have not yet attempted to combine our algorithm with multi- 
canonical methods. Hence, we expect that our cluster algorithm still suffers from 
super-critical slowing down close to the deconfinement phase transition. 

We have compared the efficiency of the Metropolis algorithm for the flux model, 
the cluster algorithm for the complex action Potts model, and the reweighted Metro- 
polis algorithm for the complex action Potts model at several points in the phase 
diagram. In figure 5 we compare the computer time histories of the Polyakov loop for 
the three algorithms mentioned above at (h, k) = (0.01, 0.5) which is in the confined 
region of the phase diagram. Obviously, the flux model Metropolis algorithm decor- 
relates a lot worse than the other two algorithms. The flux algorithm performs even 
worse when one approaches the first order transition line. The reweighted Metropo- 



20 



(a) 



J i L 



j i i i i_ 




10 20 30 40 50 60 70 80 90 



(c) 



.2 




■ n 1— 


I 1 I 

—\ ^ 





=1 1 — r 






-.2 


i 


i i i iii 


i , i 



10 20 30 40 50 60 70 80 90 

t [sec] 



Figure 5: The computer time history of the Polyakov loop $ (a) for the cluster algo- 
rithm for the Potts model using the improved estimator for (b) for the Metropolis 
algorithm for the Potts model using $ directly, and (c) for the Metropolis algorithm 
for the flux model on a 20 3 lattice at (h,K) = (0.01,0.5). The horizontal straight 
line denotes the expectation value. The Metropolis algorithm for the flux model has 
a much longer autocorrelation time than the other two algorithms. In the case of 
the Metropolis algorithm for the Potts model the complex action problem manifests 
itself by large fluctuations around the expectation value. 
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lis algorithm for the Potts model suffers from the complex action problem. What is 
plotted is the time evolution of ReE a $ x exp(ih ImQ x )]/V (exp(ih ^ x Im^^))^. 
Its statistical fluctuations are much larger than the expectation value which is a 
manifestation of the complex action problem. 



5 Universality Class of the Critical Endpoint 



In this section we present the results of our numerical simulations at the critical 
endpoint. The complex action problem turned out to be very weak in the vicinity 
of the critical endpoint (see subsection 3.4). Therefore it is possible to use a simple 
reweighted Metropolis algorithm even though that does not solve the complex action 
problem. It is usable up to sufficiently large lattices so that critical exponents can 
be extracted from a finite size scaling analysis. Of course, on even larger lattices the 
meron-cluster algorithm will eventually be superior to the reweighted Metropolis 
algorithm. But since at the endpoint E the complex action problem sets in only at 
volumes > 100 3 , simulations at E are not limited by the complex action problem 
but by the ability to simulate large lattices on today's computers. 



Figure 1 shows the phase diagram of the model defined by equation (|2.11|). For 



h = our model reduces to the standard 3-d 3-state Potts model which has been 
studied extensively in Monte Carlo simulations [pql . The model is known to have 
a weak first order phase transition. The value of the coupling k where the phase 
transition occurs (point T in fig. 1) has been determined with high precision. In 
the phase transition was found to occur at kt = 0.550565(10). Above this 



value the Z(3) symmetry is spontaneously broken, i.e. for k > kt three distinct 
deconfined phases coexist. When we switch on the parameter h, the Z(3) symmetry 
gets explicitly broken. Positive values of h favor the deconfined phase with a real 
value of (<&). Hence, the line k > kt at h = is a line of first order phase transitions 
which cannot terminate in the deconfinement transition at the point T. In fact, T 
is a triple point because two other first order transition lines emerge from it. For 
h > a line of first order transitions extends into the (h, ft)-plane and terminates in 
a critical endpoint (E in fig. 1). Negative values of h favor the two deconfined phases 
with complex values of (<&). Negative h are unphysical in the QCD interpretation 
of the Potts model because h represents exp(— (3{M — /x)) in QCD. Still, the Potts 
model at h < makes perfect sense as a statistical mechanics system (unrelated 
to QCD) and it has another first order transition line emerging from the point T. 
Interestingly, with our method the complex action problem can only be solved for 
h > since otherwise the improved estimator of eq.( |3.8|) is not necessarily positive. 
It should be noted that for h < also the flux model suffers from the complex action 
problem. 

The line of first order phase transitions Kt(h) is determined by the condition 



22 



that the free energy densities of the confined and deconfined phases are equal, 
i.e. f c {h,K t (h)) = fd(h, K t (h)). Close to the point T = (0, Kt) the free energy- 
density of the confined phase is given by 

f c {h,K) = / C)T + e c>T {K- K t ), (5.1) 

where = / c (0, k t ) and e c ^ = df c /dK,(0, «y) is the energy density of the confined 
phase at the point T. Note that to leading order f c {h, k) is independent of h because 
(<&) = in the confined phase at h = 0. On the other hand, for the deconfined phase 
one obtains 

fd{h, k) = f d , T + e d , T (K - kt) - h($) T , (5.2) 

where ($)t is the value of the Polyakov loop at the point T in the deconfined phase 
that is favored at h > 0. Using the condition / C) y = f^T for the deconfinement 
phase transition at h = 0, one finds 

($W h 

K t (h) = kt h = kt . (5-3) 

e c ,T ~ e d T r 



The Monte Carlo data of [28] and [29| imply r = 0.41(1). Our data are consistent 
with the first order transition line being a straight line. Fitting the values of K t (h) 
obtained from the infinite volume extrapolation described below yields r = 0.430(6) 
in reasonable agreement with the number from above. Similarly, one can determine 
the angle at which the third transition line leaves the point T in the direction of 
negative h. A similar argument can be applied to the Potts model with real action 
that was studied in [23j| . Also in that case the first order phase transition line is 



consistent with a straight line and again the predicted position for the transition 
line agrees with the numerical data. Interestingly, for the Potts model with both 
real and complex action, information at h = is sufficient to predict the position of 
the transition line for h > 0. This is because the transition at h = is rather weak 
and the line ends already at small values of h. If the transition would extend deep 
into the (h, «;)-plane one would expect deviations from a straight line that would be 
hard to predict based on data at h — 0. 

To determine the location of the transition line numerically we perform for given 
values of h and the volume V simulations at 3 to 5 different values of k. These sim- 
ulations are then combined with Ferrenberg Swendsen multi-histogram reweighting 
[30]. To estimate the position of the transition line we use the specific heat 



c = \ Usm - h J2 ^) 2 > - asm - h J2 ^)) 2 ) (5.4) 

\ X X / 

and determine the position of its maximum K t (h,V) for a given h and V. The 
transition point K t {h) is determined in the infinite volume limit using 

K t (h,V) = K t (h) + (5.5) 
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Figure 6: The susceptibility x along the transition line plotted as a function of h for 
five different volumes. 

where K t (h) = K t (h, V = oo). This ansatz is used sucessfully in the whole /i-range, 
i.e. for the first order region as well as in the crossover region. The values for K t (h) 
are plotted with errorbars into the phase diagram (figure 1). On the first order 
transition line they are consistent with a straight line which intersects the k axis 
exactly at the point T. The crossover line has a slight curvature. 

After we have determined the tangent to the first order transition line close to the 
endpoint we still have to find the exact location of the endpoint on that line. Also 
we want to extract critical exponents. For that purpose we consider the Polyakov 
loop susceptibility 

\ X X / 

along the transition line, i.e. at the points (K t (h), h). One could consider a "rotated 
susceptibility" instead, where an admixture of the kinetic energy term is added to 
the Polyakov loop to diagonalize the fluctuation matrix. Nevertheless, this is not 
necessary, since the "magnetic field direction" is the dominant one. We explicitly 
checked that the results do not depend on the admixture we chose, unless one comes 
close to the linear combination where the discontinuities of the Polyakov loop and 
the kinetic energy cancel. This linear combination would correspond to the kinetic 
energy term in the Ising model. Of course, if one would be interested in observables 
related to the kinetic energy of the Ising model, one would have to exactly choose 
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Figure 7: The scaling function f(x). 

the linear combination that corresponds to the kinetic energy direction. 

Close to the critical point the following scaling ansatz describes the susceptibility 
(see e.g. 0) 

X = L^f(x), x = (h-h c )L 1 'r (5.7) 

For the fit, the function f(x) is expanded in a Taylor series around x = up to 
second order, f(x) = fo + fix + f2X 2 /2. We perform two fits: 

• Six parameter fit: It results in \jv = 1.532(57), 7/^ = 2.064(70), h c = 
0.000445(18), f = 0.70(16), /1 = -3.99(95), and f 2 = 16.4(60) with % 2 /d.o.f . = 
1.33. The results for the exponents agree with the estimates for the 3d-Ising 
universality class 1/u = 1.587(2) and j/u = 1.963(3) |31[] almost within error- 
bars. 

• Four parameter fit: The critical exponents are fixed to the Ising model values. 
The result is h c = 0.000470(2), f = 0.9775(46), h = -4.567(55), and f 2 = 
16.5(17) with x 2 /d.o.f. = 1.35. The good value for \ 2 supports again the 
universality class of the 3-d Ising model. 

The susceptibility and the four parameter fit are shown in Figure 6. Figure 7 shows 
the function f(x). 
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6 Conclusions 



We have used a cluster algorithm to solve the notorious complex action problem 
in the Potts model approximation to QCD with heavy quarks at large chemical 
potential. We use a simple analytically constructed improved estimator that gives 
an exponential reduction in the required statistics. Since the improved estimator 
is real and positive, importance sampling techniques that fail for complex actions 
then become applicable. This makes it possible to study the whole h > parameter 
range of the Potts model, not just the h = axis. (Recall that h corresponds to 
exp (/?(// — M)) in QCD in the limit M, \x — > oo at any given fj, — M). 

We compared our cluster algorithm with a flux model reformulation, and with 
the reweighted Monte Carlo algorithm. We found that the cluster algorithm was 
more efficient than using the flux model reformulation. In the large volume limit 
the cluster algorithm will always be superior to Monte Carlo reweighting. However, 
at very small h it is sufficient to use reweighting techniques to obtain physically 
relevant results. This turned out to be the case for the endpoint of the first order 
line, which occurs at a very small h because the 3-d 3-state Potts model phase 
transition is rather weak at h — 0. We therefore used reweighted Monte Carlo to 
locate the first-order line and its endpoint. However, we emphasize that as computer 
power rises, and the maximum attainable volume becomes bigger, the meron-cluster 
algorithm will eventually become superior at any h > 0. 

We also calculated quark-quark, quark-antiquark, and antiquark-antiquark po- 
tentials, in the confined and deconfined regions of the phase diagram. We found the 
expected behavior: the background density of heavy quarks screens color fields, so 
that all potentials reach plateaux at long distances, whose values are simply related 
to the free energies of external static quarks and antiquarks. 

The algorithm that we have developed for the Potts model belongs to the class 
of meron-cluster algorithms that has recently been used to solve a large variety of 
sign and complex action problems. Of course, the ultimate goal is to construct a 
similar algorithm for QCD at non-zero chemical potential and investigate the phase 
structure of QCD at n ^ from first principles. The complex action problem in full 
QCD is more complicated than the one in the Potts model. So far, meron-cluster 
algorithms have led to solutions of fermion sign problems as well as complex action 
problems in bosonic theories, but have not yet solved complex action problems in 
theories with fermions. We believe that this may ultimately become possible when 
one uses the D-theory formulation of QCD. 
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